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We solve the relativistic Riemann problem in viscous matter using the relativistic Boltzmann 
equation and the relativistic causal dissipative fluid-dynamical approach of Israel and Stewart. Com- 
parisons between these two approaches clarify and point out the regime of validity of second-order 
fluid dynamics in relativistic shock phenomena. The transition from ideal to viscous shocks is 
demonstrated by varying the shear viscosity to entropy density ratio r]/s. We also find that a good 
agreement between these two approaches requires a Knudsen number Kn < 1/2. 
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I. INTRODUCTION 

One of the strongest arguments for studying a fully rel- 
ativistic formulation of fluid dynamics is the large value 
of the elliptic flow coefficient V2 measured in ultrarel- 
ativistic heavy-ion collisions at the Relativistic Heavy 
Ion Collider (RHIC) at Brookhaven National Labora- 
tory (BNL) [ij]. These measurements indicate that a 
new phase of matter, the quark-gluon plasma, is created, 
which behaves like an almost perfect fluid. An indica- 
tion for this is the good agreement between experimental 
data and theoretical calculations using ideal fluid dynam- 
ics, where the viscous transport coefficients vanish; see, 
e.g. Refs., Hii. 

In nature, the viscous transport coefficients cannot 
vanish but must have lower bounds [sUt^. Small viscous 
corrections are always needed in fluid dynamics to make 
the flow laminar and stable. It was later confirmed by 
calculations within viscous fluid dynamics [sj and micro- 
scopic transport theory that the shear viscosity has to 
be sufficiently small in order to keep the agreement with 
the V2 data. These calculations suggest that the shear 
viscosity to entropy density ratio rj/ s < 0.4. 

In heavy-ion collisions at ultrarelativistic energies the 
system expands very fast and gradients in the matter are 
very large. It is still an open question to what extent fluid 
dynamics is applicable in describing the dynamics of such 
a system. Fluid dynamics is an effective theory that de- 
scribes the macroscopic evolution of the system close to 
thermal equilibrium. Its applicability requires that either 
the viscosity or the gradients are small, or both. On the 
other hand, microscopic transport models can be used for 
systems which are also strongly out of thermal equilib- 
rium. Therefore, a comparison between the microscopic 
approach and fluid dynamics can provide the limits and 
accuracy of the fluid-dynamical description. 

In this work we will compare solutions of the relativis- 
tic Boltzmann equation with those of second-order fluid 
dynamics as derived by Israel and Stewart (IS) [lol - [T2| . 



Such a comparison has been previously presented for the 
Bjorken scaling solution [13| in Ref. [ij]. The scaling 
solution provides a simple test case with (initially) arbi- 
trarily large expansion rate. It was concluded that for a 
good agreement between the fluid-dynamical and kinetic 
calculations a Knudsen number Kn < 1/2 is required. 
The Knudsen number is the ratio of the mean-free path 
of the particles and the length scale of the variation of 
macroscopic fields (such as the energy density). For the 
same scaling solution, an extension of the IS theory was 
studied and compared to the kinetic solutions in Ref. [l^ . 

The Bjorken scaling solution is, however, restrictive 
in the sense that there are no pressure gradients in the 
local rest frame (LRF) of matter. In this paper we will 
complement these earlier works by studying the relativis- 
tic Riemann problem. This was already investigated in 
Refs. 0, [13], and we continue this line of study with 
more extensive and detailed comparisons between the 
fiuid-dynamical and transport approaches. The existing 
analytic solution in the perfect-fluid limit, i.e., for rj ~ 0, 
makes it possible to check the numerical convergence of 
the approaches. To our best knowledge this type of study 
was never accomplished for relativistic shock phenomena. 

The present work is yet another step toward under- 
standing how highly-energetic jets created in the initial 
stage of an ultrarelativistic heavy-ion collision interact 
with the hot and dense medium. If the latter is a strongly 
interacting fluid, one expects the creation of shock waves 
which resemble Mach cones. However, if matter is only 
weakly interacting, these shock waves should be smeared 
out and a clean Mach-cone signal cannot be observed 

m- 

The paper is organized as follows: In Sec. [H] the dis- 
sipative fluid-dynamical theory of Israel and Stewart is 
introduced. This is followed in Sec. Illll bv a short intro- 
duction to the numerical methods we use. In Sec. llVl the 
relativistic Riemann problem and its analytic solution 
in the perfect-fluid limit are presented. We demonstrate 
that numerical calculations using both approaches repro- 
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duce the analytic solutions very accurately. In Sec. IVlwe 
show comparisons with non-zero viscosities and demon- 
strate how the agreement between the two approaches 
deteriorates with increasing viscosity. These results are 
then analyzed in terms of a relevant Knudsen number in 
order to estimate the applicability of the IS theory. As 
in the ideal case, we have found a scaling behavior of 
the relativistic Riemann problem with non-zero viscos- 
ity. This is discussed in Sec. IVIl Finally, conclusions are 



given in Sec. IVIII 

Our units are h 
diag(l, -1,-1,-1). 



1; the metric is g 



II. FROM KINETIC THEORY TO FLUID 
DYNAMICS 

A. Definitions 

In relativistic kinetic theory of simple gases matter 
is characterized by the invariant single-particle distribu- 
tion function f(x,p). The space-time evolution of f{x,p) 
caused by particle motion and collisions is given by the 
relativistic Boltzmann transport equation [Ta . [20| , 



(1) 



where C [f{x,p)] is the collision integral and p^ — {E,p) 
is the particle four-momentum. We assume that there 
are no external forces. 

Macroscopic quantities can be obtained from the mo- 
ments of the distribution function. The first moment is 
the particle four-flow. 



N^ = J dpp'^fix,p).. 



(2) 



where dp = g d'^p/[{2Tr)^E] and g is the degeneracy fac- 
tor counting internal degrees of freedom. The second 
moment defines the energy momentum tensor, 



7^"^ = J dppy'fix^p). 



(3) 



Similarly one can define higher moments of the distribu- 
tion function, such as the third moment. 



F^^-'^ = j dpp^py'^fix^p). (4) 

The entropy four-current is defined as 

S^ = j dpp^^f{x,p) [1 - ln/(x,p)] . (5) 

The assumption of the conservation of particle number 
or charge and energy-momentum conservation in individ- 
ual collisions leads to 



dpp^dj^ dpC = 0, 



d^T'"' = / dpp^p^d^f^ / dpp'^C^O 



(6) 
(7) 



These are the conservation equations of relativistic fluid 
dynamics (2ll |. A similar equation resulting from the 
Boltzmann equation for the third moment gives the bal- 
ance of fluxes, 



I j dpp^^p-p^dxf ^ Jdpp'^p'^C. 



(8) 



This will be discussed in the next section. The Boltz- 
mann H-theorem implies that the entropy production is 
positive and vanishes in equilibrium 



d^s^ > . 



(9) 



The particle four-flow and energy-momentum tensor 
can be decomposed with respect to an arbitrary time- 
like four- vector u'^, normalized as u^^u^ = 1. This is 
chosen in such a way that it can be interpreted as the 
collective four-velocity of the matter. The frame where 
= (1,0,0,0) is called the local rest frame. With the 
help of the transverse projection operator A^'^ — g^'^ — 
u^u'^ ^ the most general decomposition can be written as 



(10) 
(11) 



where n = N^Uf^ is the LRF particle density and e = 
Uf^T^'^Ui, is the LRF energy density. The spatial trace of 
the energy-momentum tensor, P — —^/S.^^T^^'^ denotes 
the isotropic pressure. This is the sum of equilibrium 
pressure and bulk viscous pressure, P = p{e, n) +11. The 
flow of particles in the LRF is 

V = Ai^N" , (12) 

while the flow of energy-momentum in the LRF is 

(13) 

The heat flow is deflned as 

qt" = W'' -hV , (14) 

where h = (e + p)/n is the enthalpy per particle. 
The shear-stress tensor, n^'^ = T^^''\ where 



A^^Af, + AlA^;, ] - ^-A^'-A^p 



T 



(15) 



is that part of T*^", that is symmetric, traceless, and 
orthogonal to the flow velocity. 

In equilibrium all dissipative terms vanish, that is, 11 = 
= = TT^'^ = 0, therefore the particle four-current 
and energy-momentum tensor take a simpler form. 



■poA" 



(16) 
(17) 



where the subscript "0" indicates local thermodynami- 
cal equilibrium. In local equilibrium, the particle four- 
current Nq , the energy-momentum tensor Tq"^ , and the 
entropy four-current, Sq, are uniquely defined for any 
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time-like four- velocity m''. The LRF frame particle den- 
sity and energy density are given by Uq = NgUf^ and 
Co = UfjTl^'^Uv, respectively. 

For a system that deviates from local thermodynamical 
equiUbrium, the definition of the flow velocity may follow 
Eckart [l^l or Landau and Lifshitz 23]. Using Eckart's 
definition of the flow velocity, 



ATM 



(18) 



the LRF flow of particles vanishes, = 0, while the 
flow of energy is given by the heat flow, W^^ = . 

In Landau's frame the flow of matter is tied to the flow 
of energy. 



(19) 



hence the flow of energy-momentum vanishes, = 0. 

For processes close to equilibrium, the above defini- 
tions and decompositions are related to each other [23 |. 
Therefore up to second order in deviations from local 
equilibrium, that is, 5u^^ ^ g^/e ^ 1, we get 



e + p 



(20) 



This means that the non-equilibrium part of the particle 
four-flow in Landau's frame is related to the heat flow in 
Eckart's frame, that is, = —q^^/h. 



B. Fluid dynamics as an approximation to kinetic 
theory 

Fluid dynamics is an effective theory for the slow, long- 
wavelength dynamics of a given system. For systems with 
well-deflned quasi-particles, fluid dynamics can be de- 
rived in terms of a power series in the Knudsen number 



Kn 



A 



mfp 



(21) 



where Amfp is the mean-free path of the particles and 
L is a macroscopic length scale over which macroscopic 
flelds such as energy density, particle density, or temper- 
ature vary. Fluid dynamics as an effective theory can be 
systematically improved by successively including higher- 
order terms in Kn |15l] . 

To zeroth order in Kn we obtain an effective theory 
that does not contain any powers of Kn, corresponding 
to the limit Kn — > 0, that is, the (unphysical) limit where 
Amfp — >■ 0. This corresponds to infinite scattering rates, 
and thus the system instantaneously assumes local ther- 
modynamical equilibrium. This is the perfect-fluid limit. 
To flrst order in Kn, we obtain the relativistic general- 
ization of Navier- Stokes theory. This effective theory is 
plagued by instabilities and acausalities [25l - l27l |. These 
problems can be circumvented by including terms of sec- 
ond order in Kn, such as in the fluid-dynamical theory of 



Israel and Stewart used in this work. The fluid-dynamical 
limit can be derived for any kind of system, that is, its 
applicability is not restricted to dilute gases, as is the 
case for the Boltzmann equation. However, since it is an 
expansion around the perfect-fluid limit, we expect its 
validity to be restricted to dynamics close to local ther- 
modynamical equilibrium. 

The coefficients of fluid dynamics as an effective the- 
ory can be computed by matching to an underlying mi- 
croscopic theory. In our case, this is the kinetic theory 
of ultrarelativistic Boltzmann particles, described by the 
Boltzmann equation with elastic binary collisions. Note, 
however, that the matching procedure is not unique [28| . 

In our case the matching procedure is as follows. We 
expand the single-particle distribution function around 
local thermodynamical equilibrium, f{x,p) — fo{x,p) + 
Sf, where Sf measures the deviation from the equilibrium 
distribution function 



fo{x,p)^{e^P''''^/\ + ay 



(22) 



Here a = — 1(+1) corresponds to bosons (fermions), 
and a = corresponds to Boltzmann particles. The 
inverse temperature is given by I3{x) = l/T{x) and 
A(x) — e'^'^^^''^^^ is the fugacity. The validity of the ex- 
pansion around /o requires that </> = Sf/fo ^ 1. In order 
to match to fluid dynamics as an expansion in powers of 
Kn, one possibility is to assume that is a series in pow- 
ers of Kn. This approach was pioneered by Hilbert, and 
by Chapman and Enskog [20l |. 

Another method was proposed by Grad ^29"] and was 
gen eralized to relativistic systems by Israel and Stewart 
[10l - [l2| . In this approach. 



{x,p) = e + e^p^ + 



(23) 



where the coefl[icients e, e^, and e^i, are the expansion 
parameters and therefore ~ 0{Kn). 

For each non-equilibrium state given by / the corre- 
sponding equilibrium state /o is defined by the Landau 
matching conditions 12], such that riQ = n and ep — e. 
These conditions together with the Gibbs equation en- 
sure that the equilibrium part of the pressure is given by 
the equation of state po(e, n) = po(eo5 "-o)- 

The coefficients of the expansion, e, e^, and e^^^, are de- 
termined using Eq. ([^ in Eqs. © and Comparing 
the outcome with Eqs. (jTU]) and ([TT|). one finds that the 
parameters of the expansion are proportional to the dissi- 
pative quantities H, V^^, g^, and tt^''. Thus the deviation 
from equilibrium is proportional to the ratio of dissipa- 
tive quantities to local equilibrium quantities, that is, 
(f) ^ H/e, g^/e, TT^'^/e. However, if the dissipative quan- 
tities are close to the Navier-Stokes values (see later), 
then these ratios are proportional to the local Knudsen 
number. 

The macroscopic equations for the evolution of dissi- 
pative quantities can be obtained from the third moment 
of the single-particle distribution function [iH. Here we 
recall the result of this laborious calculation by Israel and 
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Stewart as presented by Huovinen and Molnar [l4|. In 
this paper we consider only massless gases in which case 
bulk viscosity vanishes. Then the IS equations for the 
heat flow and the shear-stress tensor are, 



1 



Ins ■ 



(24) 



^u^'q^Du" - (Vxu^ + Dln 



Pi Pi 



(25) 



P2 P2 

where the proper-time derivative is denoted by _D = 
the gradient operator is — A^^^d^, and the 
vorticity tensor is uj'''' = ^A^'^A'^" {dpUa - daUfi). The 
coefficients ai,/3i,/32, and oi are thermodynamic func- 
tions, and a[ = [d{Pai) / d(3) ^^j. — ai. These depend 
whether we choose the Landau or Eckart frame. The 
Navier-Stokes values for the heat flow and shear-stress 
tensor are 



^ 



»,^v^(f), 

e + p \1 / 



^NS = 



(26) 
(27) 



where Kq is the heat conductivity coefficient and rj is the 
shear viscosity coefficient. The relaxation times of heat 
conductivity and shear viscosity are proportional to the 
heat conductivity and shear viscosity coefficient, respec- 
tively, that is, Tq = KqTPi and r^r = 2?7/32. 

The microscopic time scales in the IS equations are 
given by the relaxation times of dissipative quantities 
and Tq, which are of the order of the mean- free path 
between collisions. The relevant macroscopic scales can 
be estimated from the gradients of the primary fluid- 
dynamical variables. For example, these can be given in 
terms of the expansion rate Lg — 1/9, in terms of the 
energy density gradient — ^/V^eV^e/e, or in other 
ways. 

If the Knudsen number is sufficiently small, then at 
late times t > T-„,Tq, heat flow and shear viscosity will 
approach their Navier-Stokes values, that is, q^ ^ q'^g 
and TT^'^ ~ ^jvs- When this happens, the dissipative 
quantities can be estimated to be of order 1 in the Knud- 
sen number, and Eqs. ([M)) and (j25l) include contributions 
up to second order in Kji. 

In the limit Tq^r,^ — > 0, with rj and Kq constant, the 
IS equations reduce to the Navier-Stokes equations. Fur- 
thermore, in the limit when all dissipative quantities ap- 
proach zero, the IS equations reduce to the perfect-fluid 
equations. 



C. The Israel-Stewart equations for 
(1 + 1)— dimensional expansion 

For the sake of simplicity, we assume an ultrarelativis- 
tic massless Boltzmann gas with conserved particle num- 
ber. In this case, the bulk viscosity vanishes, and the 



equation of state is simply e 



p, where the speed of 



sound is Cs — ^^1/3. For gluons, the energy density as a 
function of temperature is e = 3nT, where n — X gT^ /n^ 
is the number density, with g — 16 being the gluon num- 
ber of degrees of freedom. The entropy density is given 
by s = (4 — In X)n. 

In the following we choose the Landau frame. We 
shall briefly discuss and write the IS equations in (1-1-1)- 
dimensional Cartesian coordinates. We assume that the 
system is homogeneous in the transverse directions, x 
and y, and evolves along the longitudinal direction z such 
that the velocities as well as the derivatives in both trans- 
verse directions vanish identically. Thus the four- velocity 
is = 72(1,0,0, u^) where 7^ = (1 — v^)'^^^, while the 
four-derivative is 9^ = {dt,0,0, dz). The following four- 
vector and tensor components vanish: = — and 

j^Ox ^ yOy ^ j,xy ^ j^xz ^ j.yz ^ ^ rpj^j^ ^jg^ implies 

that the heat-flow components = = and shear- 
stress tensor components tt'^^ = ir'-^y = tt^^ = tt^^ = 
ttV^ = vanish identically. 

Using the orthogonality of the heat-flow four-vector 
we obtain that g° = q^Vz- We may also deflne the 
magnitude of the heat-flow four-vector hy q = y/—q^qfj,, 
thus q^ = "fzQ- Similarly, using the orthogonality prop- 
erty of the shear-stress tensor we get tt"*' = tt'^^Vz and 
7r°^ = TT^^Wx. To satisfy the tracelessness condition we 
may choose tt^^ = tt^^ = — 7r/2 and tt^^ = j^tt. 

Therefore, the non- vanishing components of the parti- 
cle four-current and energy-momentum tensor are 



rpOz 
fJ-iXX 



717, 



q Vz 



N%z. 



h ' 

{e + VzH-Vz.. 
VziT^^'+Vz), 

P-^ = Tyy, 

VzT°^ + Vz , 



(28) 
(29) 

(30) 
(31) 
(32) 

(33) 



where the LRF effective pressure is, 

Vz=p{e,n) + ^. (34) 

The LRF particle and energy densities expressed through 
the laboratory frame quantities and the velocity are 

-1 



q Vz 



■P 



(rpOzVJ 



T"" + Vz ' 



T 



Oz 



(35) 
(36) 
(37) 
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The conservation equations are 

dtN"" + d,{v,N'')=d.. 



q n 



(38) 
(39) 

9tr"- + 9,(z;,r"^) =-a,7',. (40) 

The relaxation equations for the heat conductivity are 
calculated from Eqs. ^ and In the (1+1)- 

dimensional case the terms containing the vorticity van- 
ish; therefore the relaxation equations can be written for- 
mally as 



9l 



— 



DtT = (tTns - Tt) - /ttI - Itt2 - Itt3 ■ 



(41) 
(42) 



where the Navier-Stokes values for the heat conductivity 
and shear stress are, 



T^NS 

Ins 



KgT\ {Ts)n 2 
s J e + p 



(43) 

dtX dzX\ , . 
+ - . (44) 



The expansion rate is denoted by 9^ = dtjz + dzi'JzVz)- 
In the ultrarelativistic limit, ai = — l/(4p), /3i ~ 5/(4p), 
p2 = 3/(4p), ai = 0, and a'l = 5ai. The terms in the 
relaxation equations are given explicitly as 



91 



T 



Iq2 = -q^Vzll (dtVz + VzdzVz) 



I. 

and 



g3 - g [iz {vzdtTT + dzir) + 72 7r {vzOz + jzdtVz 



I-irl 



0, +Dlll- 



T 



10, 



^7r2 = 77 (9^7z) (dtVz + VzdzVz) 

- I (vzdtq' + dzq' - 
9 V 7^ 



(45) 
(46) 
(47) 

(48) 
(49) 
(50) 



The terms 1^3, /7r2, and 7^3 represent a coupling be- 
tween the heat-flow four-vector and shear-stress tensor. 

In this work the term /^g is neglected in most cases 
unless otherwise stated. The reason is that the agreement 
with kinetic theory is better without it, but results with 
and without this coupling term will be shown when we 
discuss viscous solutions to the Riemann problem. 



III. NUMERICAL METHODS 

A. The transport model: BAMPS 

The microscopic transport model we use is the Boltz- 
mann approach of multiparton scatterings (BAMPS) 



[30l |3l|. This numerical method solves the Boltzmann 
equation ([T|) for on-shell particles based on the stochas- 
tic interpretation of transition rates. In this study we 
consider only binary collisions with an isotropic cross sec- 
tion, that is, a cross section with an isotropic distribution 
of the collision angle. 

Simulations of the space-time evolution of particles are 
performed in a static box. Because we consider shock- 
wave propagation in one dimension as described in the 
next section, we choose the z-axis as the direction along 
which shock waves propagate. The x- and y-axes span 
the transverse plane. 

The whole box is divided into spatial cells with a vol- 
ume V^cii — Ax Ay Az. Collisions of particles in the same 
cell are simulated by Monte Carlo technique according to 
the individual collision probability within a time step At, 



At 



Ntest Vcell 



(51) 



where a is the total cross section, and froi = (pi + 
P2)'^ / {2E1E2) denotes the relative velocity of the two in- 
coming particles with four momenta pi,P2- In order to 
reduce statistical fluctuations in simulations and to en- 
sure an accurate solution of the Boltzmann equation ^ 
a testparticle method [30] is introduced: The particle 
number is artificially increased by multiplying it by the 
number of testparticles per real particle, iVtost- Thus, the 
collision probability has to be reduced by the same num- 
ber to keep the particle mean-free path independent of 
the value of Attest- 

If a collision occurs, momenta of colliding particles are 
changed according to the isotropic distribution of the col- 
lision angle. Before and after collisions particles propa- 
gate via free streaming. Collisions of particles against 
box boundaries are realized as elastic collisions off a 
wall. We have to stop simulations before the shock front 
reaches the box boundaries, otherwise it would get re- 
flected. 

The relationship between the shear viscosity 77 and the 
total cross section a is given by 77 = 4e/(15i?''') [s^ . 
where R*^ = n(ureiCT*'') = 2n(uroif)/3 is the transport 
collision rate in the case of isotropic scattering processes 
[l4| . stands for the ensemble average in the LRF. In 
IS theory, we obtain 



2 , 



(52) 



for an ultrarelativistic massless gas, where Amfp = 
l/{n{v,:cia)) is the particle mean- free path. Note that 
(wicic) is inversely proportional to Amfp and 77/s. For 
given constant Amfp or rj/s, the average (frcic) in indi- 
vidual cells is generally different, if there exists a spatial 
gradient of n and/or e. In these cases we assume that a 
does not depend on the momenta of two incoming parti- 
cles. This leads to (frcif) « a, because (wrci) ~ 1 in the 
LRF. 

In the following, we shall also need the relationshi 
between the heat conductivity and the cross section 



i mp 

m 
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l33j which, in IS theory, is 

Kq = 2/a. (53) 

B. The viscous hydro solver: vSHASTA 

In order to solve the IS equations of causal relativis- 
tic fluid dynamics we use a version of the sharp and 
smooth transport algorithm (SHASTA) [111- This nu- 
merical method is widely used in modeling relativistic 
heavy-ion collisions, and hence was extensively tested 
in the perfect-fluid approximation [sBl, [s^. We apply 
SHASTA to solve both the conservation equations and 
the relaxation equations rearranged in conservation form 
and call this numerical method vSHASTA [37| . 

This algorithm requires the Courant-Friedrichs-Lewy 
(CFL) condition Acfl = At/Az < 0.5, where At is the 
time step and Az is the cell size. In all our numeri- 
cal calculations we take Acfl = 0.4. In a first-order 
finite-difference approach this means that causal trans- 
port of matter covers only a distance XqflAz, while the 
remaining part of the matter is acausally diffused over 
a distance (1 — Acfl)Az. This purely numerical effect 
called prediffusion is partially removed by higher-order 
non-linear corrections in SHASTA. The remaining low- 
order numerical diffusion represents the so-called numer- 
ical viscosity of the algorithm. In a strict sense numerical 
viscosity does not fully correspond to real physical vis- 
cosity, which is independent of the numerical method and 
resolution. However, its presence is inevitable and at the 
same time compulsory to keep the computations stable 
and to smooth out dispersion errors. By increasing the 
numerical resolution this numerical diffusion can be al- 
ways be reduced to smaller values than the physical one. 

We also mention another rather trivial numerical arti- 
fact which is present in the numerical solutions at early 
times. The numerical solutions at early times do not 
represent the correct and accurate physical behavior, not 
even in the perfect-fluid limit. However, this will change 
in time as the solution spreads over a larger number of 
cells while the structures are resolved on a finer grid. 
Therefore, the numerical solutions approach the correct 
solution only after some amount of time. 

IV. TESTING THE PERFECT-FLUID LIMIT 

A. The relativistic Riemann problem in the 
perfect-fluid limit 

In this section we introduce the relativistic Riemann 
problem in the case of perfect fluids. The matter is as- 
sumed to be thermodynamically normal (ssj and, for the 
sake of simplicity, to be homogeneous in the transverse 
directions, so the problem becomes (f + l)-dimensional. 

In the Riemann problem we have matter in thermody- 
namical equilibrium separated by a membrane at z — 0. 



The pressures on the left {z < 0) and right (z > 0) sides 
of the membrane are po and p4, and the particle densi- 
ties are no and n4, respectively. Here we only discuss a 
special case of the Riemann problem, called the shock- 
tube problem, when the velocities on both sides of the 
membrane are zero, that is, vq — — 0. 

Removal of the membrane at time t — leads to two 
propagating waves. If po > P4j a shock wave is propagat- 
ing to the right with the velocity Wshock- Simultaneously, 
the tail of a rarefaction fan is propagating to the left 
with the speed of sound Cs into the matter with higher 
pressure. The region between these two waves includes 
a contact discontinuity, propagating to the right with W3, 
and a shock plateau, which is bounded by the contact 
discontinuity and the shock front. 

Figure [T] shows the analytic solution for the particle 
density and velocity profile for an ultrarelativistic mass- 
less gas, e = 3p. Here, regions and 4 represent the 
undisturbed matter at rest, I is the rarefaction wave, 2 
denotes the constant region between the tail of the rar- 
efaction wave and the contact discontinuity, while 3 is 
the shock plateau. The shock front is the discontinuity 
between regions 3 and 4. 



initial conditions 




FIG. 1. (Color online) Analytic solution of the Riemann prob- 
lem for a perfect fluid, (a) particle density and (b) velocity, 
as a function of the similarity variable ^ = z/t. The initial 
temperatures are To = 0.4 GeV and T4, = 0.2 GeV. 
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The solution of the relativistic Riemann problem is ob- 
tained by matching the pressure p2 and velocity V2 at the 
rarefaction tail to the pressure and velocity of the shock 
plateau and v^, that is, p2 = Ps and V2 — v^; see Ref. 
[351] for more details. 

The solution at the discontinuity is given by the 
Rankine-Hugoniot-Taub relations [2l| in the LRF of the 
shock, 



(es + P3)73W3 = (64 + P4)l4U4 , 

(es + Ps.)ul + Pa = (e4 + Pijul + Pi . 



(54) 



Quantities with a are evaluated in the LRF of the 
shock. Then, the velocities of the shock plateau and the 
shock front in the LRF of the undisturbed matter can be 
expressed in terms of thermodynamic quantities before 
and after the discontinuity, 



^shock 



jPs -Pi){e3 - 64) 

(64 +P3)(e3 +P4) 

{P4 - P3){e3 + Pi) 



(64 - 63)(e4 +P3) 



(55) 
(56) 



where = 3pi, i = 3,4, for an ultrarelativistic gas of 
massless particles. 

The solution for the ideal shock-tube problem is self- 
similar in time, that is, the solution keeps the same shape 
at all times, t > 0, without change. This is best seen if we 
plot the solutions against the similarity variable, ^ = z/t, 
as is done in Fig. [TJ 



B. Numerical convergence of BAMPS 

Before we employ BAMPS to solve the Riemann prob- 
lem near the perfect-fluid limit in the next subsection, we 
first show the convergence of BAMPS when the model 
parameters are varied. 

In BAMPS the relevant parameters that control the 
numerical accuracy are the cell size Az, the time step 
At, and the test-particle number per real particle, iVtost- 
The cell sizes in the transverse directions, Aa; and Ay, 
are not relevant, because the system is assumed to be 
homogeneous in the transverse plane. In addition. At is 
always chosen to be smaller than Az, to avoid possible 
large local variations within one time step. If one de- 
creases Az, one has to simultaneously increase iVtcst to 
ensure that each cell contains a sufficiently large number 
of test particles. Thus, using BAMPS the Boltzmann 
equation ([T|) will be exactly solved in the limit Az — ^ 
and iVtcst — >■ 00. In practice we of course use a non- 
vanishing value of Az and a finite value of iVtost- In the 
following we show how the numerical solutions converge 
when Az is decreased and iVtest is increased. 

The initial condition for this study is Tq = 0.4 GeV 
and T4 = 0.35 GeV. Figure [5] shows the pressure profile 
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initial conditions 
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testp = 90-135 
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FIG. 2. (Color online) Cell size (a) and number of test par- 
ticle (b) dependence in the BAMPS simulation at i = 3.2 
fm/c. The initial conditions are chosen as To — 0.4 GeV 
and Ti = 0.35 GeV. The mean-free path is Amfp = 0.1 fm. 
In (a) the pressure profile is shown for different cell size 
Az = 0.4,0.1,0.025 fm and constant Nteat- In (b) we use 
a fixed Az = 0.1 fm and different numbers of test particles 
(testp) per cell. 



at a time t — S.2 fm/c for a constant mean- free path 
Amfp = 0.1 fm. Results in the upper panel are obtained 
by varying Az and keeping iVtest unchanged. The num- 
ber of test particles in the cells is between 22 and 34, 
depending on the local temperature. We see that con- 
vergence is reached when Az = Amfp. Further decrease 
of Az does not lead to noticeable changes. The lower 
panel of Fig.[2]shows the results for a fixed Az and vary- 
ing number of test particles per cell. Here we do not see 
significant changes even for a small test-particle number 
in the cells. However, as we will show in the next subsec- 
tion, a small number of testparticles in cells causes large 
fluctuations in each event, and thus affects heat flow, for 
instance. 

Results from BAMPS, which will be presented in the 
rest of the paper, are obtained by setting Az < Amfp and 
taking at least 15 test particles in each cell. We note 
that for these calculations a different initial condition. 
To = 0.4 GeV and T4 = 0.2 GeV, is chosen. For this case 
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a faster convergence has been observed. 



C. Numerical solutions of the relativistic Riemann 
problem near the perfect-fluid limit 

The Riemann problem, which is analytically solvable 
in the perfect-fluid limit, is an important test case for 
both the kinetic and the fluid-dynamical calculations. In 
this section we show that both approaches can reproduce 
the analytic solution very accurately, and also discuss 
possible numerical uncertainties. 

In BAMPS we cannot exactly reach the perfect-fluid 
limit, but we can choose a very small physical viscosity 
7] = 0.001s to simulate an ideal fluid numerically. Use of 
even smaller viscosities, or equivalently larger cross sec- 
tions, would require a better resolution (smaller Az and 
larger iVtost), which is computationally very time con- 
suming. 

On the other hand, for vSHASTA we can choose 
r] = Kq = 0, which solves the relativistic Euler equations 
instead of the IS equations. However, as explained ear- 
lier, because of the approximative nature of the numer- 
ical algorithm, we always have some residual numerical 
viscosity in the calculations. 

In both calculations the initial state was chosen such 
that To = 0.4 GeV and T4, = 0.2 GeV. Figure [3] shows 
the pressure p and velocity v, while Fig. |4] shows the LRF 
particle density n and fugacity A profiles at time t — 3.2 
fm/c. 

The BAMPS results for p, v, and n agree well with 
those of the ideal fluid-dynamical solution, except for 
small deviations around the discontinuities separating 
different regions. These deviations are expected to ap- 
pear because of the small but non- vanishing physical vis- 
cosity used in BAMPS calculations, and are best seen in 
the fugacity profile. 

Nevertheless, the BAMPS results also deviate from 
vSHASTA using the same rj/s value; this is especially 
visible in the fugacity profile between the contact discon- 
tinuity and the shock front, where one observes a peak in 
the BAMPS result. We will demonstrate later that this 
deviation is caused by numerical fluctuations. 

We also expect that particles in regions around discon- 
tinuities are out of thermal equilibrium in the BAMPS 
calculations. However, this is a small fraction of the 
whole system. The rest reaches and approximately main- 
tains thermal equilibrium. In order to demonstrate this, 
we calculate the energy distribution dN/ {N dE) of par- 
ticles in the region of the shock plateau and compare it 
to the thermal one, which is obtained by using Eq. ([^^ 
with a — Q. 



ideal hydro 
vSHASTA Ti/s = 
BAMPS r\ls = 0.001 
vSHASTA r\ls = 0.001 



> 
CD 




dN _ e-^-^/^ sinh(w7£;/T)£; 
NdE ~ 2T272W 



(57) 



Results are shown in Fig. [5j where we see agreement 
within a few percent. 
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FIG. 3. (Color online) The analytic and numerical solutions 
of the relativistic Riemann problem for the pressure (a) and 
velocity (b). The initial conditions aX t = Q are To = 0.4 GeV 
and T4 = 0.2 GeV. The full lines are the analytic solutions at 
t = 3.2 fm/c. The results from BAMPS for -q/s = 0.001 are 
shown with dashed lines, the vSHASTA results in the perfect- 
fluid limit with dash-dotted lines, and those for r;/s = 0.001 
with dotted lines. 



If the system is exactly in thermal equilibrium as de- 
scribed by the ideal fluid-dynamical solution, the dissi- 
pative quantities, such as the heat flow g^, should van- 
ish. Figure shows the heat-flow profiles from both 
vSHASTA and BAMPS calculations with the same initial 
condition with 77/s = 0.001. Heat flow in vSHASTA is 
practically zero, while in BAMPS it has a small positive 
value between the rarefaction fan and the shock front. 
This deviation between BAMPS and vSHASTA resuhs 
was already noticeable in the fugacity profile, but not in 
the other quantities shown above. At the shock front we 
see a peak, similar as for the fugacity profile. 

The deviation of from zero observed in Fig. [S] (a) 
(except for the peak) seems to be a numerical artifact, 
because this disappears when a constant cross section is 
used instead of a constant 77/s value, as shown in Fig. [S] 
(b). Here, the cross section is set to be cr = 224.431 
mb, which corresponds to 77/s = 0.002 in the medium 
with the higher initial temperature. In Fig. [B] (b) we 
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FIG. 4. (Color online) Same as in Fig. [J] The upper panel 
(a) shows the particle density while the lower panel (b) shows 
the fugacity. 



see perfect agreement between BAMPS and vSHASTA 
results, especially comparing them at the small drop at 
z « 1.6 fm. The peak at the shock front becomes smaller, 
although it is larger than the peak from the vSHASTA 
calculation. 

Also, the difference in the fugacity between BAMPS 
and vSHASTA with the same rj/s value, as observed in 
the lower panel of Fig. |H disappears almost completely 
when using a constant cross section, as seen in Fig. [T] 

On the shock plateau where the flow velocity and the 
LRF particle and energy density (w, n, e) are constant, 
one expects that the cross section is also constant for a 
constant rj/s. However, in a single event, thermodynamic 
quantities fluctuate, such that one would use a smaller 
cross section (larger shear viscosity) in a cell with larger 
energy (entropy) density in order to keep rj/s constant, 
and a larger cross section in a cell with smaller energy 
density. Therefore, although the results shown in Figs. 0] 
(lower panel) and[n](a) are averaged over 1000 events, the 
deviations between BAMPS and vSHASTA results are 
likely to come from the fluctuations in single BAMPS 
events. These can be reduced by performing simulations 
with a much larger iVtcst- We have confirmed that in 
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FIG. 5. (Color online) Relative difference of energy distribu- 
tion of particles at the plateau extracted from BAMPS and a 
local thermal equilibrium distribution, Eq. (|57|l . 



this case the difference between BAMPS and vSHASTA 
solutions decreases. 



V. VISCOUS SOLUTIONS OF THE 
RELATIVISTIC RIEMANN PROBLEM 

In this section we study the relativistic Riemann prob- 
lem at different non-zero viscosities. We will show that 
for small viscosities both the fluid-dynamical and ki- 
netic approaches are in good agreement, especially at 
late times. However, with increasing viscosity this agree- 
ment fades and ultimately breaks down when the fluid- 
dynamical description leads to results which are inconsis- 
tent with kinetic calculations. Thus the main motivation 
of this study is to find the conditions of this break-down 
and then quantify the reach and limits of the dissipative 
fluid-dynamical description. 

We note that all results shown below are calculated 
using the Landau frame. In these test cases heat flow 
is small; therefore the differences between the Landau or 
Eckart frames are very small, even for large values oirj/s. 

Kinetic theory can correctly treat the Riemann prob- 
lem from the nearly perfect limit to the free-streaming 
limit. This has been previously demonstrated in 
Refs. [l^[l3 using BAMPS. Another promising method 
to investigate the Riemann problem is based on the lat- 
tice Boltzmann approach and has been recently reported 
in Ref. [s^. In contrast to kinetic theory, the validity 
of IS theory requires that the system stays close to lo- 
cal thermal equilibrium and the Knudsen number Kn is 
small during the whole evolution. 

In the special case of the Riemann problem, at early 
times of the evolution, the local Knudsen number is large 
where the density gradients are large, even if the viscosity 
is small. Also the system around the discontinuity is far 
from equilibrium. Then in this region the IS theory of 
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FIG. 6. (Color online) The heat-flow component (a) with 
fixed rj/s — 0.001 and (b) with constant cross section a — 
224.431 mb. The value ot rj/s in the constant cross section 
simulation varies from 0.002 to 0.008. 



dissipative fluid dynamics is expected to fail to describe 
the evolution correctly. 

However, because of the viscosity and heat conduc- 
tivity the gradients will be smoothed out later, hence 
providing better conditions for the IS fluid-dynamical 
description. How close the solution will be to the ki- 
netic one depends on the value of the Knudsen number 
as demonstrated later in Sec. IV Dl 

In the next subsections we will show results at fixed 
times but for different values of rj/s. However, solutions 
at time t with shear viscous coefficient rj correspond to 
solutions at time a t with shear viscous coefficient a rj, 
where a is some arbitrary constant. This scaling behavior 
is discussed later in Sec. IVll 



A. Small viscosity 

We use the same initial conditions as given in the pre- 
vious section, but now for two different values of the 
shear viscosity to entropy density ratio, rj/s = 0.01 and 
r]/s — 0.1. Figure [8] shows the pressure p and velocity v, 
Fig. [S] shows the LRF particle density n and the fugacity 
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FIG. 7. (Color online) Fugacity calculated using a constant 
cross section a — 224.431 mb. 



A, and Fig. [TU] shows the profiles of the shear pressure 
TT and the heat flux at time t — 3.2 fm/c from both 
BAMPS and vSHASTA calculations. 

In the dissipative case the characteristic structures of 
the perfect-fluid solution can still be found in the late 
stages of the evolution, since it takes a finite time for 
the structures to form; see Sec. IVII However, instead 
of a discontinuous shock front, a contact discontinuity, 
and sharp rarefaction tails, we get continuously changing 
profiles, that is, dissipation leads to the smoothing and 
broadening of these characteristic structures. 

Further differences compared to the perfect-fiuid case 
are that the head and the tail of the rarefaction fan and 
the shock front propagate faster into the undisturbed 
matter. However, for the shock wave this happens only 
until the shock plateau is formed. After that the velocity 
of the shock wave is the same as for the perfect-fiuid case. 
Similarly, the velocity of the plateau does not change 
from the perfect-fluid solution. 

For the smaller shear viscosity to entropy density ra- 
tio, rj/s = 0.01, the agreement between BAMPS and 
vSHASTA results for all macroscopic quantities is excel- 
lent within statistical fluctuations, although any definite 
conclusions regarding the heat fiuxes are hard to draw 
because of large fluctuations. The shock front from both 
calculations is also in very good agreement. However, 
a closer inspection reveals that vSHASTA gives slightly 
steeper profiles than BAMPS. 

Increasing the viscosity to rj/s — 0.1 leads to no- 
ticeable differences between the approaches. The most 
pronounced difference can be seen in the shock front: 
vSHASTA provides a too sharp profile at the right edge 
of the front, while in the BAMPS calculation the matter 
is diffused faster in the low-density region. This can also 
be seen in the rarefaction fan where the kink at the left 
edge survives. Another difference can be seen by inspect- 
ing the fugacity and shear pressure in the region where 
the contact discontinuity would be in the perfect-fluid 
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FIG. 8. (Color online) The same as Fig. [31 for rj/s = 0.01 and 
0.1. 
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FIG. 9. (Color online) The same as Fig. [SI for (a) particle 
density and (b) fugacity. 
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case. In this region the vSHASTA calculation returns 
an overall smaller fugacity, and also a sharp kink in the 
shear pressure profile. These differences will be enhanced 
for larger viscosities as seen in the next subsection. 



B. Large viscosity 

If the value of rj/s is further increased, that is, ry/s = 
0.2, we start to see much larger deviations between 
vSHASTA and BAMPS results. These are presented in 
Fig. [TT]for the pressure and velocity profiles. 

The most salient difference is seen in the pressure pro- 
file: in the vSHASTA calculation a part of the initial 
discontinuity survives near the contact discontinuity, at 
2 ~ 1.5 fm even after t = 3.2 fm/c. In contrast, this 
kind of structure is not seen in the BAMPS calculation. 
A similar shock structure was also observed in Ref. [40l |. 
called the " double-shock" phenomenon by the authors. It 
is important to note that in that case a different numeri- 
cal method, namely, the smoothed particle hydrodynam- 
ics, was used to solve the equations of dissipative fluid 
dynamics, corresponding to the simplified IS equations 
without heat conductivity. The simplified or truncated 



IS equations only take into account the relaxation term 
to describe the time evolution of dissipative quantities. 

In the fluid-dynamical calculation this discontinuity 
originates from the initial discontinuity. In the early 
stage of the evolution the effective pressure p + n and 
velocity take almost constant values near the discontinu- 
ity. The velocity and the gradient of the effective pres- 
sure are the driving forces of the expansion. Therefore, 
if they are constant, nothing happens to the structures 
in the solution and the original discontinuity disappears 
very slowly, such that parts of it are still visible in the 
later stages of the evolution. 

A similar difference between the BAMPS and 
vSHASTA results is seen at the right edge of the shock 
wave. This part of the profile is again not described cor- 
rectly by IS fluid dynamics, and the difference is already 
visible for smaller values oirj/s. The same is true also for 
the head of the rarefaction fan, although it is less visible 
for small viscosities. This kind of discontinuous behavior 
can be seen in all relevant fluid-dynamical quantities. 

In the BAMPS calculation the original discontinuity 
disappears immediately. This is because in kinetic theory 
the evolution near the very steep density gradient is well 
approximated by free streaming or diffusion of particles. 
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which smooths out all sharp structures very rapidly. Free 
streaming of particles drives the system immediately far 
out of thermal equilibrium, hence cannot be correctly 
described by second-order dissipative fluid dynamics. 

This phenomena was studied and explained in non- 
relativistic systems M. Torrilhon et aZ.[4l|. They con- 
cluded that the viscous fluid-dynamical solutions of the 
Riemann problem actually lead to discontinuous solu- 
tions. For example, in the non-relativistic 13-field equa- 
tions the system has five instead of three characteristic 
waves. Although dissipation leads to the attenuation of 
these waves and smoothing of discontinuities, this can 
only happen after a sufficiently long time. In case we 
include higher moments of the distribution function we 
will find more characteristics and therefore more discon- 
tinuities but with smaller amplitude. There is an infinite 
number of moments, which form a hierarchy of equa- 
tions. Therefore by taking into account higher moments 
the approximations to the Boltzmann equation become 
more precise, which in turn leads to a better approxima- 
tion for smooth profiles. On this account we note that 
there are more recent studies showing that with a special 
regularization technique the Grad's 13- moment method 
leads to much better results [43|. So far these methods 



have been studied only in the non-relativistic case, but 
nevertheless they point toward a solution. 



C. Heat-flow problem 

As already mentioned in Sec. Ill CI we neglected the 
term that couples the heat flux to the shear pressure from 
the heat-flow equation (|4T]) . The reason is that, if this 
term is included, the good agreement of heat flow and 
fugacity between BAMPS and vSHASTA is lost even for 
small viscosity. This is demonstrated in Fig. 1121 where 
we show the fugacity A and heat flux q^- for rj/s = 0.1 
with and without this coupling term. The profiles change 
completely, and there is no support from BAMPS for 
structures induced by the coupling. For other quantities 
this coupling term has a very small effect. 

The reason that this single term can become dominant 
is that in the viscous Riemann problem, the heat flow is 
typically one order of magnitude smaller than the shear 
pressure. Thus, the coupling term in the heat equation 
can be large when the shear pressure is large, even if it 
is formally only a second-order correction. 
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FIG. 12. (Color online) The fugacity (a) and heat flow (b) 
profiles with and without coupling to the shear pressure in 
the heat-flow equation, 1^-^ in Eq. (|4ip 



D. Global Knudsen number analysis 

In order to better quantify and measure the applicabil- 
ity of IS theory, we define the relative difference between 
the BAMPS and vSHASTA calculations as 



Se\_J_ 
e / " Az 



dz 



Se 



(58) 



where Se is the difference in energy density between the 
BAMPS and the vSHASTA calculations. We recall that 
e — 3p for a massless gas. The integral is evaluated from 
the head of the rarefaction fan to the tail of the shock 
wave; hence the constant temperature regions to the left 
and right are not included. The width of this region is 
denoted as Az. Similarly, the average macroscopic length 
scale can be estimated from the average energy density 
gradient as 



l: 



Az I e dz 



1 , eo 

Az 64 



(59) 



where eo and 64 are the initial energy densities on the left- 
and right-hand sides of the initial discontinuity. Hence, 



an average Knudsen number relevant for this study can 
be defined as 



A 



mfp 



(60) 



where Amfp is the mean-free path in the low-temperature 
region, that is, it is the largest mean-free path. This 
definition smooths the rapid changes compared to a local 
Knudsen number and makes the comparisons between 
calculations feasible. Note that the Knudsen number ([SO]) 
is similar to that introduced in Ref. [l6j . 

Since Amfp ^ rj/{Ts) and Lf. ~ t, where t is a given 
time during the evolution, we also have 



Kup ~ — — — 
s Tt 



(61) 



Therefore, for a given temperature, Kn,, stays constant 
if we scale -q/ s and t by the same factor. 
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FIG. 13. (Color online) Time evolution of the average Knud- 
sen number, Kue, for different initial temperature ratios and 
different values of r]/s. 



Figure [T2] shows the time evolution of the average 
Knudsen number for two different initial conditions and 
two different viscosities. In all cases the Knudsen num- 
ber is initially large but decreases rapidly as the system 
expands. This happens as a function of the initial tem- 
perature difference and viscosity. 

On the other hand, the relative difference between 
BAMPS and vSHASTA calculations decreases with de- 
creasing average Knudsen number. This is shown in 
Fig. [13] where we can see that for small Knudsen num- 
bers the different solutions converge to approximately 
one curve. This also means that to good approximation 
the Knudsen number alone determines the applicability 
of IS theory. We read off Fig. [T3] that the differences 
between BAMPS and vSHASTA are less than 10% for 
Kue < 1/2. 
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FIG. 14. (Color online) The relative difference of kinetic and 
fluid-dynamical calculations for different initial temperature 
ratios and different values of r;/s as a function of the average 
Knudsen number. 



VI. FORMATION OF SHOCK WAVES AND ITS 
SCALING BEHAVIOUR 

In ideal fluid dynamics shock waves are formed imme- 
diately after removing the membrane that separates mat- 
ter with different temperatures. This happens because 
the Knudsen number Kn,, vanishes at any time. For non- 
vanishing viscosity Kue is large at early times and, thus, 
the formation of shock waves occurs later, when Kn^ be- 
comes smaller, for example, in the case for a constant 77/s 
value, as demonstrated in Fig. [T31 

Figure [13 shows the pressure and velocity profiles at 
various times for rj/s = 0.1. At the early time t = 0.64 
fm/c shock waves have not yet developed. The pressure 
profile looks like that in the strong diffusive case of free- 
streaming particles. At t = 6.4 fm/c we observe a char- 
acteristic shock plateau that clearly separates the shock 
front from the rarefaction wave, as in the ideal-fluid case. 
The intermediate time t = 3.2 fm/c is the time scale at 
which the shock plateau is being formed and the maxi- 
mum of the velocity distribution v{z) reaches the value 
Vp\nt of the ideal-fluid solution. We define this time scale 
as the formation time of shocks. 

The only intrinsic length scale in the microscopic ap- 
proach is the particle mean-free path. Therefore, if we 
rescale the mean-free path by a constant factor of a, 
we expect the time scale for the evolution of matter to 
change accordingly. Since rj/s ~ Amfp, we expect that 
profiles calculated at time t for a given value of 77/s agree 
with those at a time a t for a viscosity-to-entropy density 
ratio arj/s. This is demonstrated in Fig. 1161 where we 
show the velocity profiles as a function of the similarity 
variable i^ — zjt. 

The pressure and the velocity profile as a function of 
^ arc determined by the Knudsen number Kn,, (j60p . Ac- 
cording to Eq. ((6T|) . Kue is the same for the calculation 
with 77/s = 0.1 at t = 3.2 fm/c and that with 77/s = 0.05 



D. 




FIG. 15. (Color online) The time evolution of the shock-tube 
problem for 77/s = 0.1. The initial condition is the same as in 
Fig. El 



a.i t ~ 1.6 fm/c. Therefore, the pressure p{z,t]ri/ s)/pi) 
and velocity v{z,t;r]/s) are functions only of ^ and Kue, 
that is, p{z,t]'q/ s)/pa = F{^] Kug), and similarly for the 
velocity. For decreasing Kue the plateau of the veloc- 
ity profile in Fig. [TC] will be growing and approaches the 
shape of the ideal-fluid case shown in the lower panel of 

Fig.m 

This scaling behavior holds only for initial conditions 
with a discontinuity in pressure. If the discontinuity is 
changed to a smooth transition, the non-zero width Ftr 
of the transition region introduces another length scale. 
Then, as a function of the similarity variable, this tran- 
sition region would change under a rescaling t ^ at, 
Ctr = Ftr/i — > rtr/(at), and thus cause a different gradi- 
ent in the transition region as a function of ^. Because 
of the different initial situations, evolutions in ^ for the 
same Kue are not identical. 



VII. CONCLUSIONS 

In this work we have studied the formation and evo- 
lution of relativistic shock waves in dissipative matter 
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FIG. 16. (Color online) The scaling behavior for the shock- 
tube problem. The velocity profiles are shown as a function 
of the similarity variable = z/t for ri/s = 0.1 and ri/s — 0.05 
at times t = 3.2 fm/c and t = 1.6 fm/c, respectively. 



with non-zero shear viscosity and heat conductivity by 
solving the relativistic Riemann problem. This was ac- 
complished by using both relativistic kinetic theory and 
relativistic dissipative fluid dynamics. The relativistic 
kinetic approach solves the Boltzmann equation by us- 
ing the BAMPS code The fluid-dynamical ap- 
proach is based on Israel and Stewart theory fl3\ and 
was solved numerically by the vSHASTA method for hy- 
perbolic equations. 

After extensive comparisons between the two ap- 
proaches, we found the following: When the viscosity is 
zero, both give equivalent results. It was demonstrated 
that both approaches reproduce the analytic solutions of 
the Riemann problem in the perfect-fluid limit and the 
numerical results converge when the numerical resolution 
is sufficiently high. 

Departing from the perfect-fluid limit, for cases when 
the viscosity is small, the agreement between kinetic the- 
ory and IS theory is still excellent. As the viscosity 
increases the agreement between the approaches starts 
to deteriorate. For even larger values of the tj/s ratio, 
IS theory develops discontinuities which survive even af- 
ter long times. These sub-structures are not supported 
by the kinetic simulations. They are an artifact of the 
method of moments 41| on which IS theory is based. 
However, we also argued that part of this discrepancy 



can be understood to result from the inapplicability of 
IS theory for large Knudsen numbers. 

Quantitative statements about the applicability of IS 
theory are difficult to make in the current setup, mainly 
because the early evolution is not well described by fluid 
dynamics. This also affects the late evolution of the sys- 
tem, and therefore it is difficult to clearly separate be- 
tween the artifacts from the early fluid-dynamical solu- 
tions and the effects of large Knudsen number. 

However, we showed that a quantitative analysis in 
terms of an average Knudsen number is possible in such 
a way that it gives a good measure for the applicability 
of the IS theory. In accordance with Ref. we found 
that for Kue < 1/2 the difference between kinetic theory 
and IS theory is less than ^ 10 %. 

The shear pressure profile was reasonably well de- 
scribed by the IS equations from small to moderate vis- 
cosities. Furthermore, we also found that the results are 
quite insensitive to the second-order terms in the IS equa- 
tion for the shear viscosity. However, the same is not true 
for heat flow. At very small viscosities the heat flow was 
quite well described, including all terms in the IS equa- 
tion. At larger viscosity the results are very sensitive to 
the coupling term between the heat flow and shear vis- 
cosity. In the Riemann problem studied here, the heat 
flow is an order of magnitude smaller than the shear pres- 
sure, and the coupling to shear dominates the behavior 
of the heat flow. This coupling gives a too large contribu- 
tion in the heat- flow component, which is not supported 
by the kinetic calculations. Whether the inclusion of all 
second-order terms in the IS equations will cure the heat- 
flow problem and improve the overall applicability of the 
fluid-dynamical approach will be studied in the future. 
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